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We use the coupled cluster method implemented to high orders of approximation to investigate the 
frustrated spin-^ J 1 -J 2 -J 3 antiferromagnet on the honeycomb lattice with isotropic Heisenberg in¬ 
teractions of strength Ji > 0 between nearest-neighbor pairs, J 2 > 0 between next-nearest neighbor 
pairs, and J 3 > 0 between next-next-neareast-neighbor pairs of spins. In particular, we study both 
the ground-state (GS) and lowest-lying triplet excited-state properties in the case J 3 = J 2 = nJi, 
in the window 0 < re < 1 of the frustration parameter, which includes the (tricritical) point of 
maximum classical frustration at red = |. We present GS results for the spin stiffness, ps, and the 
zero-field uniform magnetic susceptibility, x, which complement our earlier results for the GS energy 
per spin, E/N, and staggered magnetization, M, to yield a complete set of accurate low-energy pa¬ 
rameters for the model. Our results all point towards a phase diagram containing two quasiclassical 
antiferromagnetic phases, one with Neel order for re < re^, and the other with collinear striped 
order for re > recj. The results for both x and the spin gap A provide compelling evidence for a 
disordered quantum paramagnetic phase that is gapped over a considerable portion of the interme¬ 
diate region re^^ < re < re^j, especially close to the two quantum critical points at re^^ and recj. Each 
of our fully independent sets of results for the low-energy parameters is consistent with the values 
Kci — 0.45 ± 0.02 and recj = 0.60 ± 0.02, and with the transition at re^ being of continuous (and 
hence probably of the deconfined) type and that at recj being of first-order type. 

PACS numbers: 75.10.Jm, 75.30.Kz, 75.30.Cr, 75.50.Ee 


I. INTRODUCTION 

The roles played by quantum fluctuations and frustra¬ 
tion on the ordering properties of the ground-state (GS) 
phases of Heisenberg systems of interacting spins placed 
on the sites of a regular periodic lattice continue to be 
the subject of intense study. Other things being equal, 
quantum fluctuations tend to be larger for lower spatial 
dimensionality D, lower values of the spin quantum num¬ 
ber s, and lower values of the lattice coordination number 
2 . The Mermin-Wagner theorem [l| rules out the possi¬ 
bility of magnetic order in an isotropic Heisenberg system 
when D = 1 even at zero temperature (T = 0), since it 
is not possible to break a continuous symmetry in a ID 
system even at T = 0. Similarly, the Mermin-Wagner 
theorem also implies the absence of magnetic order in 
any isotropic 2D Heisenberg model at all nonzero tem¬ 
peratures (T > 0). Thus, the study of 2D spin-lattice 
systems at T = 0 occupies a very special role for the 
study of quantum phase transitions. In order to maxi¬ 
mize the role of quantum fluctuations it is then natural 
to focus special attention on spin-^ systems on the hon¬ 
eycomb lattice, which have the lowest values of both s 
and z (= 3 in this case). 

The simplest archetypal such model is perhaps then the 


pure isotropic Heisenberg model, which comprises anti¬ 
ferromagnetic Heisenberg bonds only between nearest- 
neighbor (NN) pairs on the lattice, all with identical 
exchange coupling strength Ji > 0. This model has 
been studied many times in the past, using a variety of 
theoretical tools, including, for example, the exact di- 
agonalization (ED) of small finite-sized lattices 0 , i, 
spin-wave theory (SWT) 0], Schwinger-boson mean- 
field theory (SBMFT) 0|, a linked-cluster series expan¬ 
sion (SE) around the IsiM limit 0], quantum Monte 
Garlo (QMG) simulations (7Hlll|. and the coupled cluster 
method (CCM) [H, [l^ . From many such studies there is 
now a broad consensus that the Neel order of the classical 
(s —>■ 00 ) Heisenberg antiferromagnet (HAF) on a bipar¬ 
tite lattice at T = 0 is not destroyed by quantum fluctu¬ 
ations in the s = 5 case for the honeycomb lattice, albeit 
with a much reduced value for the Neel order parameter. 
Thus, the staggered (or sublattice) magnetization M for 
the spin-i honeycomb-lattice HAF is generally agreed to 
take a value of about 54% of the classical value. The role 
of reducing the lattice coordination number, for exam¬ 
ple, from z = 4 for the square lattice to z = 3 for the 
honeycomb lattice can be seen from the corresponding 
value for M for the spin-i square-lattice HAF, which is 
generally agreed to be about 62% of the classical value. 
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Another fundamental difference between the square 
and honeycomb lattices is that, although both are bipar¬ 
tite, whereas the former is a Bravais lattice, the latter 
is not. Instead, the honeycomb lattice has two sites per 
unit cell and comprises two interlocking triangular Bra¬ 
vais lattices. Thus, the translational invariance of the 
full honeycomb lattice is broken for any type of state in 
general. This has the consequence, for example, that the 
transition from magnetic disorder at T > 0 to Neel anti¬ 
ferromagnetic (AFM) order at T = 0 is not accompanied 
by a reduction in the spatial symmetry for the honey¬ 
comb lattice, whereas it is for the square lattice. For the 
spin-i square-lattice HAF, with one site per unit cell, 
the Lieb-Schultz-Mattis-Hastings theorem Bin then 
applies. The theorem broadly asserts that a system with 
half-odd-integer spin in the unit cell cannot have a gap 
and a nnique ground state. Thus, for the square lattice, 
any gapped state must be accompanied by a symmetry 
breaking. By contrast, for a lattice, such as the honey¬ 
comb lattice, with an even number of sites per unit cell, 
the generalization by Hastings [l^ of the Lieb-Schultz- 
Mattis theorem [l^ does not apply. In the case of spin-^ 
models on the honeycomb lattice, unlike on the square 
lattice, one can in principle have a gapped magnetically 
disordered GS phase that does not break any symmetry, 
and which has only trivial topological properties. 

It is also interesting to note that several spin-i 
honeycomb-lattice systems have been realized experi¬ 
mentally. For example, recent calculations [l^ of the low¬ 
dimensional magnetic material / 3 -CU 2 V 2 O 7 have shown 
that its magnetic properties can be described by a spin- 
i anisotropic honeycomb HAF model, albeit with two 
inequivalent NN bonds arising from the anisotropic ex¬ 
change interactions. Another example is the compound 
Na 3 Cu 2 Sb 06 in which the apparently hexagonal arrange¬ 
ment of the spin-i Cu^"*" ions in the copper oxide layers 
has been taken as evidence of honeycomb-lattice mag- 
netism da. However, it has also been pointed out [T^ 
that the structural distortion of the lattice and the or¬ 
bital states of the Cu ions again conspire to make differ¬ 
ent NN AFM bonds on the lattice sufficiently inequiva¬ 
lent as to induce ID-type magnetic behavior. Another 
compound with spin-^ Cu^“'" ions arranged in a honey¬ 
comb lattice in the copper oxide layers is InCu 2 / 3 Vi /303 
MM- This is probably the only known substance de¬ 
scribed by a spin-i HAF on the honeycomb lattice with 
equivalent NN exchange interactions. Nevertheless, even 
here decisive comparison between experiment and theory 
is made difficult by the tendency of the material to disor¬ 
der structurally, due to the mixing of the magnetic spin-^ 
Cu^’*' ions with the nonmagnetic ions [l^, [ 2 l| . 

Although the Neel order in the spin-i HAF on the 
honeycomb lattice with AFM bonds on NN sites only, all 
with the same strength Ji, is stable against quantum fluc¬ 
tuations, the much reduced value of the staggered magne¬ 
tization order parameter from its classical value suggests 
that it is likely to be rather fragile against the onset of 
frustrating interactions. In recent years, therefore, it has 


become of great interest to investigate the correspond¬ 
ing model where the NN bonds with strength Ji > 0 are 
augmented by frustrating next-nearest-neighbor (NNN) 
bonds with strength J 2 > 0 , possibly also in conjunc¬ 
tion with next-next-nearest-neighbor (NNNN) bonds of 
strength J 3 . The resulting spin-i J 1 -J 2 -J 3 model on 
the honeycomb lattice, or special cases of it (e.g., when 
J 3 = 0, or J 3 = J 2 ), have been intensively investigated 
by many authors (see, e.g.. Refs. d, refer¬ 

ences cited therein). In particular, the CCM has been 
used extensively to study the GS phase structure of the 
model [ 2 ^, In these earlier studies the system 

was mainly investigated via accurate calculations of the 
ground-state energy per spin E/N, the staggered magne¬ 
tization M, and the coefficients of susceptibility against 
the formation of various forms of valence-bond crystalline 
order. In the current paper we wish to extend the work 
to calculate both the spin gap and the complete set of 
fundamental parameters that determines the low-energy 
physics of this magnetic system. 

The low-energy properties of any strongly correlated 
system with a spontaneous symmetry breaking are gov¬ 
erned by the properties and dynamics of the correspond¬ 
ing emergent massless Goldstone bosons. For such 2D 
HAFs as are studied here, these are simply the spin- 
wave (or magnon) excitations. The interactions between 
such massless Goldstone modes, the existence of which 
in this case is due to the spontaneous breaking of the 
SU(2) spin symmetry down to its U(l) subgroup, are 
strongly constrained by symmetry considerations. A con¬ 
sistent description of the physics of the ensuing low- 
energy magnons in terms of an effective theory was pio¬ 
neered by Chakravarty et al. [^ . 

After the advent of the chiral perturbation theory 
(yPT) for the (pseudo-)Goldstone pions in quantum 
chromodynamics, a systematic low-energy effective field 
theory for mag nons was quickly developed in complete 
analogy The results obtained by yPT are exact, 

order by order, in a consistent and systematic low-energy 
expansion. They are universally applicable to models 
in the same underlying symmetry class, and where the 
symmetry is broken in the same way. The correspond¬ 
ing low-energy properties of such classes of systems are 
thus determined in terms of a (small) set of low-energy 
physical parameters that enter the effective Lagrangian 
or effective action, for example. These low-energy pa¬ 
rameters are not themselves determined by the generic 
effective field theory, but depend instead on the specific 
model. 

The fundamental low-energy parameter set that com¬ 
pletely determines in this way the low-energy physics of 
magnetic systems comprises the GS energy per particle 
E/N, the average local on-site magnetization (viz., here 
the staggered magnetization) M that plays the role of 
the order parameter, the zero-field (uniform, transverse) 
magnetic susceptibility y, the spin stiffnes s Ps, an d the 
spin-wave velocity c. In previous studies [^|^,[^ of the 
spin-i Ji-J 2 ~J 3 model on the honeycomb lattice along 
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FIG. 1. (Color online) The J 1 -J 2 -J 3 honeycomb model with Ji > 0; J 2 > 0; J 3 > 0, showing (a) the bonds (Ji = - ; 

J 2 = -;J 3 = — •—) and the Neel state, and (b) triangular Bravais lattice vectors a and b and one of the three equivalent 

striped states. Sites on sublattices A and B are shown by filled black and green circles respectively, and spins on the lattice are 
represented by the (red) arrows. 


the particularly interesting line = J 2 , we have used the 
CCM implemented to high orders to give very accurate 
calculations for the quantities E/N and M in the mag¬ 
netically ordered phases, and used them to investigate in 
detail the T = 0 phase diagram of the model. 

One of our aims here is to use the CCM now to directly 
calculate also the remaining low-energy parameters x 
Ps (and hence also he = with standard AFM 

hydrodynamics, and in units for y where the gyromag- 
netic ratio pps/h = 1). Together with our earlier work, 
a knowledge of these quantities will provide a complete 
and consistent description of the model via low-energy 
yPT. Furthermore, as we shall see, the parameters x ^-nd 
Ps themselves provide further information on the T = 0 
phase structure of the model, both for the quantum criti¬ 
cal points (QCPs) at which quasiclassical ordering melts, 
and also for the paramagnetic phases with no magnetic 
order into which the system passes. Finally, we also cal¬ 
culate the spin gap for the model, using the same CCM 
methodology, in an attempt to provide even more infor¬ 
mation about the model. 

The plan of the remainder of the paper is as follows. 
In Sec. m we first describe the model itself, including its 
classical (s —>■ 00 ) counterpart. The CCM technique is 
then briefly reviewed in Sec. Ell where we concentrate on 
its key features, before presenting our detailed numerical 
results in Sec. IIVI We end in Sec. IVlwith a summary and 
conclusions. 


II. THE MODEL 

The J 1 -J 2 -J 3 model on the honeycomb lattice is de¬ 
scribed by the Hamiltonian 

H = Ji'^Si-Sj+J2 Si-Sk + Ja ^ Si • S; , (1) 
do> (d,fe» ((d,0» 

in which the index i runs over all N lattice sites, and the 
indices j, k, and I run respectively over all NN, NNN, 


and NNNN sites to i, and where in each sum each bond 
is counted once and once only. Each site i of the lattice 
carries a spin-i particle described by the spin operator 
Si = (sf, s^,sf). The lattice and the Heisenberg exchange 
bonds are illustrated in Fig. [TJ In the present paper we 
will be interested in the case where each of the three types 
of bonds is AFM in nature, i.e., J„ > 0; n = 1, 2, 3. 

The honeycomb lattice is bipartite with two triangular 
Bravais sublattices A and B. If the lattice spacing on the 
honeycomb lattice (i.e., the distance between NN sites) 
is d, then sites on sublattice A are at positions = 
ma-l- nb = '/%{m — ^n)dx + ^ndz, in terms of real-space 
basis vectors a = y/Zdx and b = id(— + 3z) for 
the honeycomb lattice, which is defined to lie in the xz 
plane, as shown in Fig. [TJ Each unit cell i at position 
vector R.; comprises a pair of sites, situated at R^ € 
A and (Ri + dz) G B. The parallelogram formed by 
a and b defines the honeycomb Wigner-Seitz unit cell. 
The Wigner-Seitz cell may itself equivalently be taken as 
being centered on a point of sixfold symmetry so that it 
is bounded by the sides of a primitive hexagon of side 
d. The first Brillouin zone is then itself also a hexagon, 
which is rotated by 90° with respect to the Wigner-Seitz 
hexagon, and which has a side of length 47r/(3-\/3d). 

The classical version of the honeycomb model of Eq. 
(HD (i .e^when s —>■ 00 ) has been studied b y s everal au¬ 
thors [22l-l^. For example, Rastelli et al. searched 
for coplanar or uniformly canted spin configurations that 
minimize the classical energy and found the former to 
be energetically favored. Generally they correspond to 
spiral configurations described by a wave vector Q, plus 
an angle 0 that describes the relative orientation of the 
two spins in the same unit cell i, both of which are now 
characterized by the same triangular Bravais lattice vec¬ 
tor Ri. The classical spins (of length s) in unit cell i are 
thus given by 

Si,r = s[cOs(Q • Ri -I- 9r)Xs + Sin(Q • R^ -I- 9r)Zs] , (2) 

where Xs and Zs are two orthogonal unit vectors that 
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define the plane of the spins in spin space, as shown in 
Fig. m The index r labels the two sites in the unit cell. 
Clearly we may choose the angles Or so that = 0 
for spins on sublattice A and 9b = 9, say, for spins on 
sublattice B. 

In our regime of interest (i.e., when J„ > 0; n = 1,2, 3) 
the classical model has a (T = 0) GS phase diagram 
comprising three distinct phases [2J, [2^ . One value of 
the spiral wave vector that minimizes the classical GS 
energy is 


Q 



■ Jl - 2 J2 ■ 

_4(J2 - J3)_ 


(3) 


together with 9 = tt. Glearly, the wave vector Q of Eq. 
© is only properly defined when 


- 1 < 


Jl - 2 J2 

4(J2 - J 3 ) 


< 1 , 


or, equivalently, when 


(4) 


3 1 11 

y<^x--] y<-x + -, ( 5 ) 

where x = J2/J1 and y = J3/J1. The region in the 
positive quadrant (i.e., x > 0, y > 0) of the xy plane that 
satisfies both inequalities of Eq. m is where the classical 
Ji“J2“ J3 model on the honeycomb lattice has the spiral 
phase described by the wave vector Q of Eq. ([3|) , and 9 = 
TT, as the stable GS phase. Along the line y = \x—j, Q = 
r = (0,0), and the phase described by Eq. ([2]) simply 
becomes (continuously) the collinear Neel AFM phase 
illustrated in Fig.[IJa). The phase transition between the 
Neel and spiral phases is thus a continuous one. Similarly, 
along the line y = 52 ;+Q = 27r/(•\/3d)i, and the phase 
described by Eq. @ becomes (continuously) the collinear 
striped AEM phase illustrated in Eig. [IJb). The phase 
transition between the striped and spiral phases is thus 
also a continuous one. The above two phase boundaries 
meet at the point (a;,y) = ( 5 , 5 ) which is the classical 
tricritical point. There is a first-order phase transition 
between the two collinear states (i.e., the Neel and striped 
states) along the boundary line x = ^, for y > i. In 
summary, when Ji > 0 and x > 0, y > 0, the classical Ji- 
J2-J3 model on the honeycomb lattice has three stable 
GS phases (at T = 0): a Neel AFM phase for y > |x — j, 
^ < X < ^ and y > 0, 0 < X < |;a striped AFM 
phase for y > ^x + j, x > i; and a spiral phase for 
0<y<|x — 4, i<a;<4 and 0 < y < ^x -I- 4, x > 4. 

We note that the spiral and the striped states described 
by Eq. ([3]) and 9 = tt have two other similar states rotated 
by in the honeycomb plane. When x —>■ 00 for a fixed 
finite value of y (i.e., when the J 2 bond dominates), the 
spiral pitch angle (^ = cos“^[|(Ji — 2 J 2 )/(J 2 — J 3 )] —>■ | 7 r. 
In this limiting case the classical model thus becomes 
two HAFs on disconnected interpenetrating triangular 
lattices with the classical 3-sublattice Neel ordering of 
NN spins oriented at angle | 7 r to each other on each sub¬ 
lattice. In this limiting case the wave vector Q of Eq. ([3]) 


becomes one of the six corners of the first Brillouin zone, 
= 47r/(3-\/3(i)x. Glearly, there are only two distinct 
such corner vectors, and these describe the two inequiva¬ 
lent 3-sublattice Neel orderings for a classical triangular 
HAE. 

We also note that when the spiral pitch angle (p takes 
a value in the range |7r < </) < tt the wave vector Q 
of Eq. Q lies outside the first Brillouin zone. It can 
equivalently be mapped back into the first Brillouin zone 
in this case, when Q moves continuously from a corner 
= —2tt/{ 3'/3d)x + 2TT/{3d)z of the Brillouin zone 
along one of its edges to the midpoint = 2'K/{3d)z. 
The striped state shown in Eig. [IJb) may thus be equiv¬ 
alently described by the ordering wave vector Q = 

(with the relative angle between sublattices A and B 
given by 0 = tt). The other two striped states are thus 
given by the wave vectors of the two other inequivalent 
midpoints of edges of the first Brillouin zone, = 

TT/{'/3d)x + Tr/{3d)z and = —7r/(-\/3(i)x -I- n/{3d)z 
(and in both of these cases with 0 = 0). 

Although the spin ordering of Eq. usually suffices 
to find all classical GS configurations , there is an as¬ 
sumption that the GS order is either unique (up to the 
trivial degeneracy associated with a global spin rotation) 
or exhibits, at most, a discrete degeneracy. However, ex¬ 
ceptions can arise for special sets {Q} of wave vectors 
[ 23 , [MI . These include the case when Q is equal to one 
half or one quarter of a reciprocal lattice vector G. This 
is precisely the case for the striped states where the wave 
vectors Q = mW, i = 1,2,3 are just one half of cor¬ 
responding reciprocal lattice vectors. As explained by 
Eouet et al. [ 2 ^ the GS ordering in this case spans a 2D 
manifold of non-planar GS configurations, all degenerate 
in energy. 

It is well known that classical models that exhibit such 
an infinitely degenerate family (IDF) of GS phases in 
some region of the T = 0 phase space often lead to the 
emergence of novel quantum phases in the correspond¬ 
ing quantum-mechanical model. A typical scenario then 
finds that quantum fluctuations lift this (accidental) GS 
degeneracy, either completely or in part by the well- 
known order by disorder mechanism [5ll - l53l| , so that just 
one or several members of the classical IDF are favored. 
Indeed, thermal or quantum fluctuations do select the 
collinear striped states out of the 2D IDF manifold, ac¬ 
cording to Ref. [ 2 ^ , wherein it is also explicitly shown in 
an ED study of the finite-lattice spectra that the degen¬ 
eracy is lifted in favor of the collinear ordering. 

In the extreme s = 5 quantum limit one may also ex¬ 
pect that quantum fluctuations might be strong enough 
to destroy any quasiclassical magnetic long-range order 
(LRO) completely in some region of the T = 0 GS phase 
space. The goal of finding any such novel quantum phases 
with no magnetic LRO, and delimiting their region of sta¬ 
bility in the T = 0 phase space, has provided the impetus 
for much recent work i [ 22 I - I 4 I I on the spin-4 J 1 -J 2 -J 3 
model on the honeycomb lattice. A particularly chal¬ 
lenging, yet potentially fruitful, regime is to consider the 
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case J 3 = J 2 = kJi since it includes the point of max¬ 
imum classical frustration, viz., the tricritical point at 
K = Kci = Henceforth, therefore, we restrict ourselves 
to this regime. 

In our earlier paper [ 2 ^ we presented results for the 
GS energy and magnetic order parameter of the spin- 
i J1-J2-J3 model on the honeycomb lattice along the 
line J 3 = J 2 = kJi, with Ji > 0 and k > 0 , using the 
CCM implemented in high orders. We found that the 
first-order transition in the classical (s —^ 00 ) model at 
K = Kci = ^ between the Neel and collinear striped AFM 
phases is split into two transitions for the s = 5 model. 
The Neel phase was found to survive for k < Kci ~ 0.47, 
and the striped phase for n > Kc^ ~ 0.60. In the re¬ 
gion Kci < K < Kc 2 between the two quasiclassical phases 
we found a paramagnetic phase with no discernible mag¬ 
netic order. By further calculating within the same CCM 
methodology, the susceptibilities of the two AFM phases 
against the formation of a state with plaquette valence- 
bond crystalline (PVBC) order, we concluded that the 
intermediate state was one with PVBC order. The ev¬ 
idence from those calculations was that the quantum 
phase transition (QPT) at Kc 2 appears to be a first-order 
one, while that at Kc^ is of continuous type. Since the 
Neel and PVBC phases break different symmetries, we 
concluded that the quantum transition at between 
these two phases is of the deconfined type . 

Our aim in the present paper is to shed further light on 
the model by calculating other physical properties within 
the same CCM methodology as used previously. Firstly, 
in order to gain more evidence about the nature of the 
intermediate phase we calculate the spin gap, i.e., the en¬ 
ergy gap between the ground state and the lowest-lying 
(magnon) triplet excitation. Secondly, as discussed pre¬ 
viously, we also calculate for both quasiclassical phases 
the spin stiffness, ps, and the zero-field (uniform, trans¬ 
verse) magnetic susceptibility, y, in order both to provide 
a complete set of low-energy parameters for both phases 
of the model and to use these parameters to provide more 
numerical evidence for the two QPTs at Kc^ and Kc 2 • 

The spin stiffness (or helicity modulus) of a spin-lattice 
system is a measure of the energy required to rotate 
the order parameter of a magnetically ordered thermody¬ 
namic system by an (infinitesimal) angle 9 per unit length 
in a given direction. Thus, if E{6) is the GS energy as 
a function of the imposed twist and N is the number of 
lattice sites, we have 


m 

N 


E{e = 0) 
N 




( 6 ) 


Note that 6 has the dimensions of an inverse length. In 
the thermodynamic limit {N —^ 00 ) a nonzero (positive) 
value of Pa implies that the system has magnetic long- 
range order (LRO), while the magnetic LRO melts at 
the point where ps 0. Clearly, for the Neel AFM state 
of Fig.HKa), for which the ordering wave vector Q = F = 
( 0 , 0 ), the quantity ps is independent of the direction of 
the applied twist, whereas for the particular striped AFM 




FIG. 2. (Color online) The twisted reference states for the 
calculation of the spin stiffness coefficient, ps, for the Ji- 
J 2 -J 3 honeycomb model, showing the twist applied in the x 
direction to the (a) Neel and (b) striped states. The spins on 
lattice sites • are represented by the (red) arrows. 


state shown in Fig. [TJb), for which Q = 2tt/{\/ 2,d)x, the 
physically relevant direction in which to apply the twist is 
the X direction. Figure [2] thus illustrates the two twisted 
AFM states (i.e., the Neel and striped states) used in 
our calculations for the spin stiffness coefficient. The 
definition of Eq. © easily leads to the corresponding 
values of ps for the classical (s — 00 ) J 1 -J 2 -T 3 model on 
the honeycomb lattice, 

= J(^i-6J2+4J3)d2s2, (7) 

and 

pfr^ = l{-Ji-2J2+U3)dh\ ( 8 ) 

The lines along which the spin stiffness coefficients given 
by Eqs. o and dH) vanish are, as expected, just the 
classical phase boundaries for the two AFM states with 
the spiral state. 

It is also interesting to note that in the limiting case 
J 3 —>■ 00 , for fixed values of Ji and J 2 , the Ji~J 2 ~ 
J 3 model reduces to four decoupled honeycomb-lattice 
HAFs, each with NN coupling J 3 and lattice spacing 2d. 
Thus, the GS energy in the case {Ji = 0, J 2 = 0, J 3 = 1} 
should be equal to that of the case when {Ji = 1, J 2 = 
0, J 3 = 0}. Similarly, the spin stiffness in the former limit 
should equal four times that in the latter limit, due to the 
doubling of the lattice size of each of the four decoupled 
honeycomb sublattices. Just, as this result holds for the 
classical model (s —?► 00 ), from Eq. 0 , so it should also 
hold for general values of the spin quantum number s. 

In order to calculate the zero-field magnetic suscepti¬ 
bility, X, we now place our system in an external trans¬ 
verse magnetic field h. For the two collinear AFM states 
shown in Fig. [U both with spins aligned along the Xs 
direction, we apply the held in the Zs direction, h = hzs- 
The Hamiltonian H = H{h = 0) of Eq. (IT|) thus becomes 

H{h) = H{h = 0) + hJ2st, (9) 

I 
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FIG. 3. (Color online) The canted reference states for the 
calculation of the zero-field magnetic susceptibility, y, for the 
Ji-J 2 “J 3 honeycomb model. The external magnetic field is 
applied in the Za direction to the (a) Neel and (b) striped 
states, shown in Fig. [T] The spins on lattice sites • are rep¬ 
resented by the (red) arrows. 


in units where the gyromagnetic ratio = 1. In the 

presence of the field, the spins will now cant at an angle 
a to the Xs axis with respect to their zero-field configura¬ 
tions, as shown in Fig.[3]for the two classical AFM states 
illustrated in FiglTJ The classical (s —>■ oo) value of a 
for each of the states is easily calculated by minimizing 
the classical energy, E = E{h), corresponding to Eq. dH]) 
with respect to a. The uniform magnetic susceptibility 
is then defined, as usual, by 

'W-vS' <“> 

and its zero-field limit as y = x(0). The corresponding 
analog of Eq. is thus. 




( 11 ) 


These relations readily yield the values of x for the two 
collinear AEM states of the classical (s —>■ oo) Ji~J 2 -J 3 
model on the honeycomb lattice. 


..Neel 

Xcl 


1 


6( Ji + J3) 


( 12 ) 


and 


striped 

Xcl 


1 

2 ( Ji -h 4J2 -h 3J3) ’ 


(13) 


with both parameters independent of s in the classical 
case. The two values of x become equal (but nonzero, 
note) along the line J 2 = (independent of J 3 ), which 
is the classical phase boundary between the two AEM 
states. 

We note that our definitions for both ps and x are per 
unit site, as is more usual for a discrete lattice descrip¬ 
tion. On the other hand, in continuum field-theoretic 
terms, such as in descriptions using yPT, it is more nat¬ 
ural to define corresponding quantities, ps and y, per 
unit area. Since the honeycomb lattice has 4/(3\/3d^) 
sites per unit area, ps = j\/3d^ps and y = jy/3d‘^x- 


III. THE COUPLED CLUSTER METHOD 

The CCM is one of only a very few universal methods 
of ab initio quantum many-body theory (QMBT) that are 
capable of systematic improvement in well-defined hierar¬ 
chies of approximations. It is nowadays one of the most 
pervasive, most powerful, and most accurate at attain¬ 
able levels of computational implementation, of all fully 
microscopic formulations of QMBT. It has been applied 
with considerable success to a wide range of physical sys¬ 
tems (ssl - ls^ , ranging from the electron gas to atoms and 
molecules of interest in quantum chemistry, from nuclear 
matter and finite atomic nuclei to strongly interacting 
quantum field theories, and from models in quantum op¬ 
tics, quantum electronics, and solid state optoelectronics 
to diverse condensed matter systems. The CCM has thus 
yielded accurate numerical results for a very wide range 
of both finite and extended systems defined either in a 
spatial continuum or on a regular discrete lattice. Of 
particular relevance to the present paper is the fact that 
the method has already been applied with demonstra¬ 
ble success to a large number of spin-lattice problems 
of interest in quantum magnetism (and see, e.g.. Refs. 
[ 1 ^ [ 1 ^ [ 2 ^, 135^^ l60l - [70l| and references cited therein). 

Two characteristic and rather unique features of the 
CCM that are immediately worthy of mention are: (i) 
its ability to deal with infinite [N ^ 00 ) systems from 
the outset, hence obviating the need for any finite-size 
scaling, which is required in most competing methods; 
and (ii) the fact that it exactly preserves the important 
Hellmann-Eeynman theorem at all levels of approximate 
implementation. The method is implemented in practice, 
as we explain below, at various levels of approximation 
in a well-defined truncation hierarchy, with each level 
specihed by a truncation index n = 1, 2, 3, • • •. The only 
approximation made is then to extrapolate the values 
obtained for the physical observables at the nth level to 
the n —?► 00 limit where the CCM becomes exact. 

The CCM has been described in detail in earlier papers 
(and see, e.g., Refs. [H, ISTl - f^ and ref¬ 

erences cited therein), and hence we briefly outline only 
its most salient features here. Every implementation of 
the method begins with the choice of a suitable normal¬ 
ized reference (or model) state, with respect to which the 
quantum correlations present in the exact CS phase un¬ 
der study can then be incorporated at the next stage. 
Eor this study suitable choices of the model state |$) will 
be the two quasiclassical AEM states (viz., the Neel and 
collinear striped states) shown in Figs.[TJa) andlljb). 

The exact GS ket- and bra-state wave functions, which 
satisfy the respective Schrbdinger equations, 

H\'i>) = E\4>) ] {^\H = E{^\, (14) 

are chosen to have normalization conditions such that 
= (4)|'I') = (4>|<I>) = 1. These exact states are 
then parametrized in terms of the respective model state 
as 

|vl') = e®|<f>); (4-1 = (<f>|5e-^. 


(15) 
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with the exponential parametrization being a key char¬ 
acteristic feature of the CCM. The two correlation oper¬ 
ators are then themselves formally decomposed as 

S = Y,SiC+; S=1 + J2 SiCY, (16) 

7#0 1^0 

where we dehne Cq = 1 to be the identity operator, 
and where the set index I denotes a complete set of 
single-particle configurations for all of the particles. In 
our present spin-lattice application it defines a specific 
multispin-flip configuration with respect to the model 
state 1 $), such that is the corresponding wave 

function for this configuration of spins. The model state 
|d>) thus acts as s fiducial (or cyclic) vector or, in other 
words, as a generalized vacuum state, with respect to 
the complete set of mutually commuting creation op¬ 
erators {C/}, and which hence satisfy the conditions 
($1(7/ = 0 = (771$), VI 7 Oj where the destruction 
operators = {Cfy. 

It is very convenient for spin-lattice problems to con¬ 
sider each lattice site i as totally equivalent to all others, 
whatever the choice of model state |$), and one sim¬ 
ple way to do this is to make a passive rotation of each 
spin so that in its own local spin-coordinate frame it 
points, say, in the downward (i.e., negative Zg) direction 
as in the spin coordinate frame shown in Fig. [I] Such 
choices of local spin coordinates clearly do not affect the 
basic SU(2) spin commutation relations. However, all 
independent-spin product model states now take the uni¬ 
versal form 1$) = I 4 , 4 , 4 , ••• 1). Thus, (7/ can simply 
be expressed as a product of single-spin raising opera¬ 
tors, s/ = -h isl, such that (7/ = s/s/^ ’ ’ ’ = 

1, 2, • • • , 2sN. For the present study, where we consider 
s = |, each site index included in the corresponding set 
index I = {ki,k 2 , - ■ ■ ,kn] n = 1, 2, • • • , 2sN} may ap¬ 
pear no more than once. Once the local spin coordinates 
have been chosen for the particular model state |$), the 
Hamiltonian H simply needs to be re-expressed in terms 
of them. 

In principle, what remains is then to calculate the 
CCM correlation coefficients {Si, Si}. This is achieved 
by minimization of the GS energy expectation functional, 

77 = i7[5/,5/] = ($|5e-^i7e^|$), (17) 

with respect to each of the coefficients {Si, Si ;V/ 0}. 

From Eqs. (HH) and (HI), variation with respect to Si 
yields the coupled set of non-linear equations, 

($|C7e"^i7e^|$) =0, V/7 0 , (18) 

for the coefficients {5/}. Similarly, variation of Eq. (flTll 
with respect to Si yields the corresponding set of linear 
equations 

($|5e"^[i7,(7/]e'^|$) = 0, V/ 70 , (19) 

for the coefficients {5/}, once Eq. (fTSl) has been solved 
for the coefficients {5/}. The GS energy eigenvalue E, 


which is the value of H at the minimum, is then simply 

E = ($|e“'^iJe^|$) = ($|iJe^|$). ( 20 ) 

Equation dUD may then be written in the equivalent 
form, 

($|5(e-^i7e^-£;)(7/|$) =0, V/7 0 , (21) 

of a set of generalized linear eigenvalue equations for the 
coefficients {Si}. 

Having built in the correlations into our CCM 
parametrization of the GS wave function |$), it now suf¬ 
fices to apply a linear excitation operator X® to |$) to 
parametrize the excited-state wave function |$e) Si'S 

|$e) = 7^^e®|$), = XfCf . ( 22 ) 

//o 

By suitably combining the GS Schrbdinger equation (fHl) 
with its excited-state counterpart, 

H\'\!,) = Ee\^), (23) 

and realizing that the operators 77® and S commute, by 
construction, we readily deduce the equation, 

e-^[H, 7s:®]e®|$) = Ae7s:®|$), (24) 

where = E^ — E is the excitation energy. By taking 
the overlap of Eq. (IMl) with the state ($|(77 we find 

($|(77e-^[Ff,77®]e'5|$) = AeT-/, V/ 70 , (25) 

when the states labelled by the indices {/} are, as usual, 
orthonormalized, ($|(77(7/|$) = S{I,J). The general¬ 
ized eigenvalue equation (1551) is then solved for Ag. 

In the present case the configurations chosen in the 
expansion of Eq. (l22l) for the excitation operator A® are 
restricted to those that change the z component of the 
total spin, S^, by one. The value we thus obtain for Ag 
is then the spin-triplet gap, which we henceforth denote 
as A. 

Up to this point no approximations have been made. 
Nevertheless, the equations (1151) for the coefficients {Si} 
are intrinsically nonlinear and one may wonder if trun¬ 
cations are needed in the evaluation of the exponential 
functions. We note that these appear, however, only in 
the combination of the similarity transform e~^Ele^ of 
the Hamiltonian, which may be expanded in terms of the 
well-known nested commutator sum. Another key fea¬ 
ture of the CCM is that this otherwise infinite sum now 
actually terminates exactly with the double commutator 
term. This is due to the basic SU(2) commutation rela¬ 
tions and because all of the terms in Eq. comprising 
S commute with one another and are simple products 
of spin-raising operators. All terms in the expansion of 
e~^He^ are thus linked, and the Goldstone linked cluster 
theorem is exactly preserved, even if the expansion of Eq. 
m is truncated in any way, thereby guaranteeing size- 
extensivity at any such level of truncation and our ability 


(29) 


to work from the outset in the thermodynamic {N —^ oo) 
limit. Similar considerations also apply to Eqs. (HU) and 

(EH). 

Thus, the only approximation made in practice for the 
GS calculation is to restrict the set of indices {/} retained 
in the expansions of Eq. (fTHl) for the operators {S', 5}. 
Here we utilize the well-tested localized (lattice-animal- 
based subsysteml LSUBn scheme used in our earlier work 
on this model and in many other studies too. 

The LSUBn scheme is defined so that at the nth level of 
approximation we retain all multispin-flip configurations 
{/} in Eq. (I16II that are defined over n or fewer contigu¬ 
ous lattice sites. Such cluster configurations are said to 
be contiguous in this sense if every site in the cluster is 
NN to at least one other. The number, Nf = Nf{n), of 
such distinct fundamental configurations may be reduced 
by utilizing the space- and point-group symmetries of the 
model, together with any conservation laws that pertain 
to both the Hamiltonian and the specific model state be¬ 
ing used. Even so, the number Nf(n) increases rapidly 
as the truncation index n is increased, and the need even¬ 
tually arises to use massive parallelization together with 
supercomputing resources for the highest-order calcula¬ 
tions [^,(7ll- 

Eor the excited-state calculation of the spin gap A the 
choice of clusters for the excitation operator A® which 
are retained in Eq. is different from that of the GS, 
since we restrict ourselves now to those that change the 
2 ; component, S^, of total spin by one unit. Neverthe¬ 
less, we use the same LSUBn scheme for both the GS 
and excited-state calculations, thereby ensuring compa¬ 
rable accuracy for both. The number Nf = Nf{n) at 
a given level n of truncation is appreciably higher for 
the excited-state than for the corresponding GS calcu¬ 
lation. Nevertheless the corresponding GGM equations 
have been solved for the present model using both quasi- 
classical (Neel and striped) AEM states as model states 
in LSUBn approximations with n < 12. 

Clearly, the CCM LSUBn approximations become ex¬ 
act, by construction, in the n —>■ 00 limit. There exist 
well-tested, accurate extrapolation rules for the GS quan¬ 
tities E/N and M, as we have described and used in our 
earlier paper for this model. Similarly, for the spin 
gap, we use the extrapolation scheme [TOI. IT^. 

A(n) = do + dirT^ + d 2 n~^ , (26) 

to obtain the extrapolated value A = A (00) = do 
from the CCM LSUBn approximations, A(n). Similar 
schemes have also been successfully used previously for 
both the spin stiffness ps [H, , 

Ps{n) = So + sin~^ + S 2 n~'^ , (27) 

and the zero-field magnetic susceptibility, y [1^ , 

x{n) = xo + xin~^ + X2n~'^ . (28) 

In the latter case, as a check on the validity and accuracy 
of the scheme, we also utilize the completely unbiased 


scheme, 

X(n) = xo +xin ~'', 

in which the leading exponent v is itself a fitting parame¬ 
ter, along with xq and xi. Finally, since the lowest-order 
LSUBn approximants (particularly that with n = 2) are 
less likely to conform well to these extrapolation rules 
than those with higher values of n, and also since the 
hexagon is the fundamental structural element of the 
honeycomb lattice, we prefer to use LSUBn data with 
n > 6 , whenever practicable, to perform each of the ex¬ 
trapolations in practice. 

IV. RESULTS 

In our earlier work on the spin-i J1-J2-J3 model on 
the honeycomb lattice along the line J 3 = J 2 (= nJi) 
in phase space [ 2 ^ we employed the CCM and com¬ 
puted LSUBn results for both the GS energy per spin, 
E/N, and magnetic order parameter, M, with values of 
the truncation index n < 12, using both the Neel and 
striped collinear AFM states as CCM model states. Al¬ 
though the number of fundamental configurations, Nf, 
at a given LSUBn level of approximation is greater for 
the triplet excited state than for the ground state for the 
corresponding CCM calculations based on both model 
states, we are still able to calculate the spin gap A at 
LSUBn levels with n < 12, even with the increased com¬ 
putational difficulty. We are thus able to achieve compa¬ 
rable accuracy for both the ground and excited states. 

We show in Fig. [4{a) our “raw” LSUBn results for 
the spin gap A with n = {6,8,10,12}, based on both 
the Neel and striped collinear states used separately as 
the CCM model state. We also show the correspond¬ 
ing (LSUBoo) extrapolated values, as obtained using the 
scheme of Eq. ((26l) . For the Neel results our method of 
solution is as follows. We first solve for the pure Heisen¬ 
berg AFM in the case k = 0, where we find a stable 
physical solution to the CCM equations at each LSUBn 
level of truncation. For a given value of the truncation 
index n the corresponding LSUBn solution is tracked as 
the frustration parameter k is increased incrementally, 
just to the point where this stable solution terminates, 
as shown in Fig. [41(a). These termination points for the 
excited-state CCM equations are completely analogous to 
those also found for the GS CCM equations, which have 
been well described and documented previously (see, e.g., 
Refs. [H, m, [g^, [g^). They are direct manifestations of 
the respective QCP in the physical system under study, 
at which the corresponding form of magnetic order in 
the associated model state melts. As is usually the case, 
we find that each LSUBn Neel solution with a fixed (fi¬ 
nite) value of n terminates at a higher value of k than 
the corresponding actual critical value, Kcj , which is just 
the LSUBoo limiting value. The outcome is that we may 
thereby consider a range for the parameter k that is ap¬ 
preciably beyond the (seemingly continuous) transition 
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FIG. 4. (Color online) The spin gap A versus the frustration parameter k = J 2 /J 1 , for the spin-i Ji~J 2 -J 3 model on the 
honeycomb lattice (with Ji = 1, J 3 = J 2 > 0). (a) CCM results based on the Neel and the striped AFM states as CCM model 
states. The LSUBn results are shown with n = {6,8,10,12}. The extrapolated LSUBoo results using Eq. (I26II are shown, 
together with the error bars associated with the fit. (b) ED results for lattices with N = 24 and N = 32 spins. 


at Kci from the Neel phase to the quantum paramagnetic 
phase. 

On the other side of the phase diagram we similarly 
track each LSUBn solution based on the collinear striped 
state as CCM model state from high values of k down to 
some respective lower transition point, corresponding to 
the actual transition at Kc 2 ■ Again, for each value of the 
transition index n, we may enter into the region below the 
(seemingly first-order) transition at Kc^ into the quantum 
paramagnetic phase. Of course, if the transition at Kc 2 is 
indeed of first-order type, as seems likely from all prior 
available evidence, one might possibly query the validity 
of our CCM results based on the striped collinear state 
in the region k < Kc 2 ■ 

In Fig. |4}a) for the LSUBoo extrapolated values of A 
based on our LSUBn results with n = (6,8,10,12}, we 
also show the error bars associated with the assumed 
scheme of Eq. (1^51) . It is clear that the fit on the Neel 
side of the phase diagram is markedly better than that on 
the striped collinear side. Very interestingly, very similar 
behavior was also observed in a recent CCM study of the 
spin gap in the spin-i J 1 -J 2 model on the square lattice 
[70| , where possible reasons for the difference in accuracy 
of the fits for small and large values of n were discussed. 
Despite this difference in the quality of the extrapolations 
for the Neel and striped phases, the results in Fig. |4}a) 
are clearly compatible with A being zero in the ranges 
K < Kci and K > Kc 2 for the two quasiclassical GS phases 
with magnetic LRO, where the spin gap must be zero. 

Conversely, we see from Fig. |4}a) that A > 0 for a 
considerable range of values in the range Kd < k < Kc 2 
on both the Neel {k > Kci) and striped {k < Kcj) sides of 
the region, over which the LSUBn solutions with n < 12 
persist, before their respective termination points. The 
values of k at which the LSUBoo curves for A become 
nonzero are also completely compatible with the values 
for Kci and Kc 2 determined by us previously [ 2 ^ , at which 


the magnetic order parameter vanishes, M —>■ 0. Our re¬ 
sults are incompatible with the entire intermediate phase 
for Kci < K < Kc 2 being a gapless spin liquid. By contrast 
they provide supporting evidence to our earlier findings 
[ 2 ^ that the quantum paramagnetic phase in this regime 
has PVBC order. On the other hand, our results for A 
would not, by themselves, rule out a gapped spin liquid 
in the intermediate regime. 

It is interesting to note that our results for A being 
nonzero over (at least part of) the intermediate regime 
for the present honeycomb model, are in contrast to the 
equivalent CCM results [t^ for the corresponding inter¬ 
mediate quantum paramagnetic regime in the spin-i Ji- 
J 2 model on the square lattice. For the latter case the 
extrapolated values of A over the entire parameter re¬ 
gion accessible were found to be zero (or very close to 
zero). For the spin-i J 1 -J 2 model on the square lattice 
the Neel order is found to melt at a comparable value of 
the frustration parameter, k= J 2 /J 1 , at Kc^ ~ 0.45, with 
a paramagnetic state in the region Kci < k < Kc 2 ~ 0.59. 
In this case the CCM based on the Neel state as model 
state can access the region Kr, < k < 0.49 in LSUBn ap¬ 
proximations with n < 12 [70|, and in this region around 
the seemingly continuous transition at Kc^ the extrapo¬ 
lated CCM value of A is zero within computational accu¬ 
racy. A very recent, accurate, density-matrix renormal¬ 
ization group calculation also finds good evidence for 
a near-critical state in the region Kc^ < ^ ^ 0.5, with a 
very small gap for the hnite-sized systems studied, and 
also for a gapped PVBC state in the remainder of the 
paramagnetic regime, 0.5 < k < Kc 2 - 

Finally, on the issue of the spin gap, we also present 
corresponding results in Fig. |3}b) for the current model 
obtained by ED, for lattices containing N = 24, 32 spins. 
Clearly, neither lattice is large enough to show clearly 
the A —>■ 00 behavior of A = 0 over both quasiclassi¬ 
cal regimes (i.e., Neel and striped collinear). Neverthe- 
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less, there is evidence (interestingly, more marked for the 
smaller lattice) of a gap opening up around k rs 0.6 (i.e., 
at Kca) as K is decreased from above. The steep decay 
around this value (seen most clearly in the TV = 32 data) 
is due to a level crossing in the triplet state. While the 
ED data are clearly indicating the first-order transition 
at K = Kc 2 , they are, unsurprisingly, quite smooth around 
the continuous transition at k = Kci ■ Thus, while the ED 
results certainly do not contradict our CCM results, by 
themselves they are certainly much less predictive than 
the corresponding CCM results. Certainly, they do not, 
by themselves, contradict the appearance of a gap in the 
intermediate regime /tci < k < Kc 2 with a sizable peak 
value, A rs 0.4, as shown by the CCM results. 

We turn now to our results for the low-energy pa¬ 
rameters, ps and X- la view of the lower symmetries 
of the twisted and canted CCM model states used re¬ 
spectively for these parameters, as illustrated in Figs. [2] 
and 131 the numbers Nf of fundamental configurations 
at a given LSUBn level of approximation are greater 
than those for the GS parameters E/N and M and 
for the spin gap A. For example, for the calculation 
of A, Nf Ri 10^ (2 X 10®) for LSUBn approximations 
based on the Neel model state with n = 10(12). Cor¬ 
responding values based on the striped model state are 
Nf Ri 2.5 X 10^ (5 X 10®) with n = 10 (12). By contrast, 
for the twisted model states used for the calculation of 
Ps, as shown in Fig. [2l TV/ r; 3.5 x 10® at the LSUBIO 
level of approximation, and it becomes computationally 
infeasible to perform LSUBn calculation for ps at higher 
truncation levels, n > 12 . 

Figure [5] shows our “raw” CCM LSUBn results with 
n = { 6 , 8 , 10 } for ps, together with the corresponding 
LSUBoo values using the extrapolation scheme of Eq. 
(l27l) together with this data set. For comparison pur¬ 
poses, we also show in Fig. [5] the corresponding clas¬ 
sical values obtained from Eqs. © and ([8|) by putting 
.h = J 2 = i^Ji and s = i, viz., p^^il{Ji(P) = ^( 1 - 2 k), 
and = ^(~1 + 2 k). We see clearly 

that the extrapolated values for ps for both (the Neel 
and striped collinear) quasiclassical phases lie somewhat 
lower than their classical counterparts for all values of 
the frustration parameter k. Based on LSUBn extrap¬ 
olations of Ps with n = { 6 , 8 , 10 }, the critical values at 
which Ps —>• 0 on the Neel and striped sides of the phase 
diagram, respectively, are k^ ~ 0.433 and Kc 2 ^ 0.621. 
These may be compared with the corresponding values 
in our earlier paper [2^ , based on LSUBn extrapolations 
of the corresponding points where the magnetic order pa¬ 
rameter TU —>• 0, of Kci R! 0.448 and ~ 0.601 based 
on the same data set n = { 6 , 8 , 10 }, and the presumably 
even more accurate values k^ R 0.466 and Kc^ R 0.601 
based on the larger data set n = { 6 , 8 , 10 , 12 } available 
in this case. Clearly, the totally independent estimates 
of the two QCPs from the places where TU —>■ 0 and 
Ps —> 0 are in good agreement with one another, within 
very small errors arising from the extrapolations. 

In order to estimate the accuracy of our results for ps 
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FIG. 5. (Color online) CCM results for the spin stiffness ps 
(in units of J\d^) versus the frustration parameter k = J 2 /J 1 , 
for the spin -1 model on the honeycomb lattice (with 

Ji = 1, J 3 = J 2 > 0). We show results based on the Neel 
AFM and the striped AFM states as CCM model states. The 
LSUBn results are shown with n = { 6 , 8 ,10}, together with 
the extrapolated LSUBoo results using Eq. ([27ll with this data 
set. The classical results from Eqs. © and ((S]) are also shown 
for the value s = |. 


independently, we may compare with those of others for 
the special case of a pure honeycomb-lattice HAF with 
NN interactions only, i.e., when k = 0. Our LSUBoo 
extrapolated value using the LSUBn data set with n = 
{6,8,10} is Ps(k = 0) r 0.1324Ji(T^. To our knowledge 
the best available alternative result for Ps(k = 0 ) comes 
from a first-principles QMC method using a highly effi¬ 
cient loop-cluster algorithm [^, [TJ . The most accurate 
value quoted by Jiang m is Ps{k = 0) = 0.1012(2)Ji, 
equivalent to Ps(k = 0) = 0.1315(3)Jid^, which is in 
excellent agreement with our own result. 

Figure [5] shows our corresponding results for the zero- 
field magnetic susceptibility %. For the striped collinear 
state used as the CCM model state, the number TV/ 
of fundamental configurations at the LSUBIO level is 
TV/ R 3.5 X 10®, and once again, just as for ps it is 
computationally infeasible to perform LSUBn calcula¬ 
tions for this phase with n > 12. By contrast, for the 
Neel state used as our CCM model state, TV/ r 6 x lO'^ 
at the LSUBIO level and TV/ r 1.1 x 10® at the LSUB12 
level. In order to be consistent with our treatment of 
the two quasiclassical phases, we show the “raw” LSUBn 
results in Fig. [5] with n = {6,8,10} and corresponding 
LSUBoo extrapolations based on this data set. However, 
on the Neel side we have also performed an LSUB12 cal¬ 
culation for the unfrustrated limiting case k = 0. Once 
again, for comparison purposes, we also show in Fig.IHlthe 
corresponding classical values obtained from Eqs. ca) 
and (fTSl) by putting J 3 = J 2 = kJi and s = ^, viz., 
= l/[ 6 Ji(l + k)], and = l/[2Ji(l + 7k)]. 

Once again we see that quantum fluctuations act to 
reduce the value of x considerably from its classical value 








11 



1C 


FIG. 6. (Color online) COM results for the zero-field magnetic 
susceptibility (in units of versus the frustration param¬ 
eter K = J 2 IJ 1 , for the spin-i J 1 -J 2 -J 3 model on the hon¬ 
eycomb lattice (with Ji = 1, J 3 = J 2 > 0). We show results 
based on the Neel and the striped AFM states as COM model 
states. The LSUBn results are shown with n = {6, 8,10}, to¬ 
gether with the extrapolated LSUBoo(l) and LSUBoo(2) re¬ 
sults using Eqs. (I28II and (1291) . respectively, with this data set. 
The classical results from Eqs. m and (O are also shown 
for the value s = ^. 


in both quasiclassical AFM phases. Figure | 6 ] also shows 
that the two extrapolation schemes of Eqs. (1^ and (1^ 
give results in very close agreement with one another on 
both sides of the phase diagram, except in very narrow 
regions close to the two QCPs at k = and k = Kc 2 - 
Perhaps the most important feature of Fig. [5] is that, 
unlike in the classical case where Xci takes the (nonzero) 
value i at the phase transition point k^i = in the 
quantum s = ^ case x now vanishes at the two QCPs, 
which is a very clear indication of a spin gap opening up 
at these points [zl, [75| . 

We show in Fig. [5] extrapolated LSUBoo results based 
on the schemes of both Eqs. (E5)) and (1^ . The two sets 
of results shown are seen to be in excellent agreement 
with each other except in very narrow regimes close to 
the two QCPs at Kci and Kcz- Based on the scheme 
of Eq. (I28L the two critical values obtained from ex¬ 
trapolating the LSUBn values with n = {6,8,10} are 
Kci ~ 0.469 and Kc 2 ~ 0.583. The corresponding val¬ 
ues from using scheme of Eq. (Ei) are Kci ~ 0.430 and 
Kc 2 ~ 0.609. While both sets of results agree quite well 
with the corresponding results cited above of Kc^ ~ 0.448 
and Kc 2 ~ 0.601 from our earlier work on extrapola¬ 
tions for the magnetic order parameter of corresponding 
LSUBn results with n = ( 6 , 8 ,10}, there is clearly more 
uncertainty in the critical values obtained from the zero- 
field magnetic susceptibility. 

As an aside here, we remark that the shapes (e.g., the 
slopes) of the respective extrapolated (LSUBoo) curves 
for ps and x in Figs. [5] and [5] appear to vary significantly 
near where they vanish at the QCP Kci, and also, but to a 


lesser extent, at the QCP Kcj . Clearly, in turn, this would 
imply different values for the critical indices for the two 
low-energy parameters at in particular. However, 
we remain cautious about putting too much weight on 
such an interpretation, since it is precisely in the regions 
very close to the QCPs where the extrapolations become 
most demanding. While we are rather confident of the 
calculated values for Kc^ and Kc 2 themselves, within er¬ 
rors we can estimate (and also see the further discussion 
below in Sec. 13), the actual detailed shapes of the extrap¬ 
olated curves (including, e.g., their slopes) at the QCPs 
is more open to doubt and the errors are more difficult to 
quantify. A close comparison of the raw LSUBn families 
of curves and their respective LSUBoo extrapolations in 
Figs. [S] and H] shows that the slopes of the raw and ex¬ 
trapolated curves very near to a supposedly continuous 
transition like that at Kci can differ appreciably. This 
difference is less marked at a first-order transition like 
that at Kc 2 , although still present to a lesser degree. For 
these reasons we are reluctant, with our present method¬ 
ology and accuracy, to make any claims to be able to 
calculate, with any quantitative degree of accuracy, the 
critical indices for the vanishing of ps and x at the two 
QCPs. Any such interpretation drawn from our results 
about the critical indices being different for ps and x at 
Kci especially, should be regarded as suggestive at best, 
in our opinion. 

It is also of interest to compare our results for x for 
the limiting case k = 0 of a pure honeycomb HAF with 
NN interactions only. Our extrapolated LSUBoo results 
based on Eq. (l28l) are x(«^ = 0) = 0.0847(4)/Ji, us¬ 
ing LSUBn data points n = {6,8,10,12}, where the 
quoted errors is purely that associated with the fit. Cor¬ 
responding extrapolations based on LSUBn results with 
n = { 6 , 8 , 10 } and n = { 8 , 10 , 12 } are, respectively, 
x(k = 0) a; 0.0845/Ji and x(k = 0) a; 0.0837/Ji. 
Other estimates are x('^ = 0) ~ 0.0756(10)/Ji from 
a linked-cluster SE analysis Q; x(^ = 0) ~ 0.1667/Ji 
and x(«^ = 0) ~ 0.0456/Ji from SWT at leading (clas¬ 
sical) order and next-to-leading order [viz., with 0 (l/s) 
corrections included], respectively and x(^ = 0) = 
0.0666/Ji from SBMFT Q. The most accurate esti¬ 
mates for this unfrustrated limit undoubtedly come from 
QMC calculations. For example. Low 0 ' used a con¬ 
tinuous Euclidean time QMC algorithm to find a value 
x(k = 0) = 0.05188(8)/Ji. Presumably, in the QMC 
calculations of the magnetic susceptibility the value x 
obtained is an average over all directions of the applied 
field, since the ground state is calculated for a finite 
system, i.e., with no breaking of the rotational sym¬ 
metry. By contrast, in the CCM calculations we start 
with a symmetry-broken state and apply the field per¬ 
pendicular to the axis of the order parameter, to ob¬ 
tain x± (= x) directly. Thus, x = 5 (X|| + 2X-l), where 
X|| (= 0) is the parallel component of the magnetic sus¬ 
ceptibility. Hence, we have x = |Xi and the QMC result 
of Low is equivalent to x(^ = 0) = 0.0778(1)/Ji. An¬ 
other QMC estimate, using a loop-cluster algorithm by 
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FIG. 7. (Color online) COM results for the three renormal¬ 
ization constants corresponding to the spin stiffness pa, the 
zero-field magnetic susceptibility x, and the spin-wave veloc¬ 
ity, c, of the spin-i Ji~J 2 ~J 3 model on the honeycomb lattice 
with Ji > 0, in the case J 3 = J 2 = kJi. Results for Zc 
are unreliable in the regions shown by thin lines, close to the 
critical points (and see text for details). 



Jiang [HI, can also be quoted. While Jiang does not 
quote a result for x directly, he does provide the result 
hc{K = 0) = 1.2905(8) Jid. This may be combined with 
his result for ps cited above to calculate x = ps/{hcY, 
yielding the value x(«^ = 0) = 0.0789(3)/Ji. Clearly, our 
own best CCM estimate, x(«^ = 0) = 0.084(2)/Ji, lies 
slightly higher than those two QMC estimates. 

It is also convenient to express the low-energy param¬ 
eters in terms of a multiplicative renormalization con¬ 
stant Z with respect to the corresponding classical result. 
Thus, we define Zp^ and Z^ as 

Zp^ = Ps/ Ps;cl j Zy_ = x/Xci • (30) 

Since the spin-wave velocity is given by he = \fp7Jx^ its 
corresponding renormalization constant is 

= \l^Pa /-^X • (31) 

Our extrapolated LSUBc» results for the three renormal¬ 
ization constants of Eqs. (1301) and (IHTl) are shown in Fig. 
[3 where we employ the LSUBoo(l) results for x obtained 
from the extrapolation scheme of Eq. (1^51) . 

We should point out that the seeming vanishing of the 
Goldstone spin-wave velocity at both critical points 
and Kc 2 , as shown in Fig. [3 is entirely an artefact of 
our calculational scheme. Thus, Zc has been calculated 
here indirectly via Eq. dMD using our (separately) ex¬ 
trapolated results for ps and x- In reality both of these 
parameters should go to zero at the same points (at least 
for X at a transition at which a spin gap opens). How¬ 
ever, as we have noted, our calculated results for ps and 
X are completely independent of one another, and as a 
consequence the critical values at which both parameters 
vanish are not constrained to be the same. This is both a 


strength and a weakness of the methodology. The main 
advantage is that it provides an inbuilt error estimate 
for the accuracy of our calculated values of and Kc 2 ■ 
However, it is also clearly a disadvantage when we wish to 
take ratios, as in Eq. dSU, in regimes close to the QCPs. 
Thus, if Ps and x go to zero at slightly different values of 
K, it is guaranteed that the spin-wave velocity calculated 
in terms of them will, quite artificially, approach either 
zero or infinity near the actual QCP. We have thus indi¬ 
cated in Fig. [3 by thinner portions of the curves for Zc, 
those regions near the two QCPs where the results for Zc 
are correspondingly unreliable. 

Finally, we collect in Table |T] our CCM results for the 
full set of low-energy parameters for the unfrustrated 
(k = 0) limiting case of a pure spin-l HAF on the hon¬ 
eycomb lattice with NN interactions only (of strength 
Ji >0). We compare our results there with values ob¬ 
tained from a linked-cluster SE analysis @ , and two dif¬ 
ferent QMC analyses which are expected to be 

very accurate in this case, where the infamous minus-sign 
problem is absent. There is no reason to expect that the 
evident accuracy of our results in this limit will be any 
lower over the entire range of values of k accessible to us, 
by strong contrast with QMC techniques, which degrade 
significantly in the presence of frustration. 


V. SUMMARY AND CONCLUSIONS 

In this paper we have continued our prior investigation 
[ 2 ^ of the spin-1 HAF model on the honey¬ 

comb lattice, along the line J3 = J2 = kJi, with Ji > 0, 
that includes the point of maximum classical frustration 
at K = Kci = 5 . Just as in our prior work we have used 
the CCM based on both the Neel and collinear striped 
AFM states as reference states, with respect to which 
we have included quantum fluctuations in the fully con¬ 
sistent LSUBn truncation scheme. We have carried out 
calculations to high orders in the truncation index n (typ¬ 
ically with n < 10, but in some cases also with n < 12 ). 
The only approximation made has been to extrapolate in 
the truncation index n to the exact {n —> c») LSUBoo 
limit. 

We have now calculated a complete set of low-energy 
parameters {E/N,M,ps,x, and c) for the model, from 
which we obtain independent pieces of evidence that pro¬ 
vide a clear and consistent description of its T = 0 quan¬ 
tum phase diagram. In particular, we find compelling 
evidence that the single phase transition in the classical 
(s —>■ 00) version of the model at KcI = ^ is split in the 
s = 1 model into two quantum phase transitions at Kci 
and Kc 2 ■ In our prior work these two QCPs were iden¬ 
tified independently by the points at which both the mag¬ 
netic order parameter M and the inverse of the suscep¬ 
tibility coefficient against the formation of a state with 
PVBC order, 1/xp, vanish. To those estimates we now 
add two more, based on the points at which both pg and 
X vanish. The collected results are displayed together in 
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TABLE I. The low-energy parameters for the spin-| HAF on the honeycomb lattice with lattice spacing d and with exchange 
interactions between NN pairs only, all with equal strength Ji > 0. The classical values are compared with our CCM results, 
and those using the alternative techniques of linked-cluster SE and two different QMC algorithms Em. 


Parameter Classical 

value 

CCM 

SE 

(Ref. [ 6 ]) 

QMC 
(Ref. [1^) 

QMC 
(Ref. [1^) 

A/(AJi) 

-0.375 

-0.54466(2) 

-0.5443(3) 

-0.54455(20) 


M 

0.5 

0.2714(10) 

0.266(9) 

0.2681(8) 

0.26882(3) 

Ps/{Jid^) 

0.1875 

0.1324(5) 



0.1315(3) 

JiX 

0.1667 

0.084(2) 

0.0756(10) 

0.07782(12) 


hc/{Jid) 

1.0607 

1.255(15) 



1.2905(8) 


Tableini in which the results for M and l/xp are obtained 
[ 2 ^ from LSUBcxd extrapolation based on LSUBn data 
points with n = { 6 , 8 , 10 , 12 }, while those for ps and x 
are obtained from data points with n = ( 6 , 8 ,10}. All of 
our results are fully consistent with values Kci = 0.45(2) 
and Kc 2 = 0.60(2). 

Based on the shape of the curves for M and 1/xp 
as functions of k we suggested in our earlier work that 
the QCP at between the state with Neel order (for 
K < Kci) and the paramagnetic intermediate state was 
of continuous type, while that at Kc 2 between the state 
with collinear striped order (for k > KC 2 ) and the inter¬ 
mediate state is of first-order type. The new evidence, 
based on ps and Xi shown in Figs.[S]and|n]respectively, is 
completely consistent with this interpretation, as indeed 
is also the evidence based on the curves for the spin gap 
A shown in Figs.|3}a) and|3}b). 

Since at a QCP the quantum fluctuations present in 
a system are sufficiently strong to make the system in¬ 
finitely susceptibility to multiple forms of order, the van¬ 
ishing of 1 /xp at Kci and Kc 2 cannot be taken as strong 
evidence that PVBC order is present over the entire in¬ 
termediate regime, Kci < k < Kc 2 - However, just as we 


TABLE 11. Values of the two quantum critical points and 
Kc 2 of the spin-i J 1 -J 2 -J 3 HAF model on the honeycomb lat¬ 
tice, along the line J 3 = J 2 = kJi with (Ji > 0), as obtained 
from the vanishing points of the magnetic order parameter 
M, the spin stiffness Ps, the uniform magnetic susceptibility, 
and the inverse of the susceptibility coefficient of the system 
against the formation of PVBC order 1/Xpi all evaluated by 
the extrapolation of CCM LSUBn results. 


Parameter (—>■ 0) 


^C2 

M (see Ref. [2^) 

0.466 

0.601 

ps (this work) 

0.433 

0.621 

X (this work): LSUBoo(l) “ 

0.469 

0.583 

X (this work): LSUBoo(2) ** 

0.430 

0.609 

1/Xp (see Ref. 12911 

0.473 

0.586 


“Using the extrapolation scheme of Eq. (I28II . 
*'Using the extrapolation scheme of Eq. (I29II . 


could in this work find values of A into this regime in 
regions around both QCPs, so in our earlier work [ 2 ^ 
could we calculate Xp in similar regions. The shape of 
the curves for 1/xp as a function of k was consistent 
with 1/xp vanishing in those regions. It was on this ev¬ 
idence that we made the tentative conclusion that the 
entire region Kr, < k < Kr„ contained a quantum phase 
with PVBC order. 

The results of the present paper provide strong support 
for such an interpretation. Firstly, the vanishing of x at 
Kci and Kc 2 provides compelling evidence for a gapped 
state opening up at these points. Secondly, our results 
displayed in Fig. |4] provide positive and conclusive evi¬ 
dence of a gapped state over a considerable range of the 
intermediate regime, which is accessible using both qua- 
siclassical AFM states as CCM model states. The values 
at which A becomes nonzero are also wholly consistent 
with the values for and Kc 2 cited above, as obtained 
from the GS low-energy parameters. 

We note that Goldstone’s theorem implies that any 
state that breaks spin-rotational symmetry must have a 
vanishing gap. Thus, the non-vanishing of the triplet gap 
A in the intermediate phase completely rules out the pos¬ 
sibility of any type of magnetic order being present in this 
regime, not only the Neel and striped forms considered 
explicitly here. Similarly, the fact that A 0 in the in¬ 
termediate phase also precludes other more exotic forms 
of order that break SU(2) symmetry. Examples include 
spin-nematic states, which break SU(2) symmetry while 
preserving translational and time-reversal symmetries. 

Finally, we note that it might also be interesting in fu¬ 
ture work to calculate the singlet excitation gap within 
the disordered regime, in order to compare it with the 
triplet gap calculated here. While the singlet gap can 
certainly also be calculated within the CCM framework, 
the calculations are more challenging, since the excited 
state now lies in the same sector as the ground state, 
viz., with = 0 , where ^ possible 

motivation for doing so would be to compare with the 
corresponding results for the spin-i J 1 -J 2 HAF on the 
square lattice, the phase diagram for which is qualita¬ 
tively similar to that for the spin-i J 1 -J 2 -J 3 Heisenberg 
model on the honeycomb lattice, in the case J 3 = J 2 
considered here. 
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For example, a recent highly accurate density-matrix 
renormalization group (DMRG) simulation of the spin-i 
J 1 -J 2 Heisenberg model on the square lattice [t^ found 
that the singlet gap remains consistently below the triplet 
gap over the intermediate disordered regime in this case. 
The authors took this finding as an indication of the for¬ 
mation of short-range singlets in the intermediate phase. 
In turn, this would be consistent with the intermediate 
phase being a spin liquid or one with only weak valence- 
bond crystalline (VBC) order. By contrast, a phase with 
stronger VBC order would be expected to have a triplon 
excitation as the lowest-energy excited state, since such 
a state corresponds to the breaking of only one singlet 
bond, compared to a singlet excitation that requires the 
breaking of two singlet bonds. If our conclusion for the 
present case of a spin-i J 1 -J 2 -J 3 Heisenberg model on 
the honeycomb lattice (with J 3 = J 2 ), that the interme¬ 
diate phase has PVBC order, is correct, we would then 
expect the singlet excitation gap to lie higher than the 
triplet gap, by contrast with the above DMRG findings 
for the spin-i HAF on the square lattice. 

In conclusion, the present paper revisits a model to 
which the CCM had previously been applied, with the 
joint aims (and outcomes) to improve the conceptual 
framework and to yield new physics, particularly with 


regard to the nature of the intermediate phase. Thus, 
the calculation of a complete set of low-energy parame¬ 
ters for the model within a single, unified, and consistent 
theoretical GGM framework, has not only given more de¬ 
tailed information about each of the ordered magnetic 
phases and more accurate values for their boundaries 
with the intermediate (non-classical) disordered phase, 
but has also opened the possibility for a full (quantita¬ 
tive) xPT treatment of the model. Similarly, the GGM 
calculation of the triplet gap A has now definitively ruled 
out the intermediate phase from having any form of or¬ 
der that breaks SU(2) symmetry, including such exotic 
states as spin nematics. 
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